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ABSTRACT 



This thesis analyzes the dynamic stability of positively buoyant submersibles. Six 
degree -of-freedom equations of motion are used to compute steady state behavior with 
motion restricted to the vertical plane. Steady state solutions are analyzed for various 
conditions of buoyancy including changes in (I) the amount of excess buoyancy, (2) the 
location of the center of buoyancy, (3) the location of the center of gravity, as well as (4) 
the deflection of bow and stern planes. The equations of motion are then linearized around 
these steady state solutions to predict dynamic response in the vertical plane. The stability 
of each solution is then determined by eigen value analysis. The study then expands the 
analysis to include all six degrees of freedom (i.e., include stability analysis in the 
horizontal plane). Finally, numerical integration methods arc used to verify the results. 
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I. INTRODUCTION 



A. PROBLEM STATEMENT 

Controlling emergency ascent situations on submersible vehicles such as dive plane 
jam recovery is ol' concern to the U.S. Navy. In order to control such situations, one must 
first be able to predict the dynamic response of positively buoyant submersibles. 

Dynamic response equations of motion descibe the manuevering characteristics of 
submersible vehicles for six degrees of freedom. These equations assume constant 
coefficients for hydrodynamic forces and moments approximated by zero frequency added 
mass and damping terms plus the quadratic terms for drag forces. The constant coefficients 
vary for each vehicle and are dependent on such things as vehicle body shape, location and 
magnitude of vehicle weight, location and magnitude of vehicle buoyancy, position of bow 
and stern planes, position of rudder, vehicle speed, vehicle mass characteristics, vehicle 
hydrodynamic coefficients, propeller rpm and control surface inputs. This thesis uses the 
equations of motion and hydrodynamic coefficients for a submerged Mark IX Swimmer 
Delivery Vehicle (SDV) developed by Smith, Crane, and Summey [Reference 1:11-16,21- 
31 ] to forecast the dynamic stability of a submersible in a positively buoyant condition. 

This study begins by using the six equations of motion to compute the steady state 
behavior of a submersible vehicle with motion restricted to the vertical plane. The steady 
state solutions in the vertical plane are calculated for various conditions of buoyancy 
including changes in the amount of excess buoyancy, the location of the center of 
buoyancy, the location of the center of gravity, as well as the dellcction of bow and stern 
planes. The SDV's equations of motion are then linearized around these steady state 
solutions to predict dynamic response motion in the vertical plane for the various conditions 



of buoyancy. Several solutions are computed and the stability of each solution is 
determined by eigen value analysis. The thesis then expands the analysis to include all six 
degrees of freedom (i.c. include stability analysis in the horizontal plane). Finally, 
numerical integration methods are used to verify the results. 

B. EQUATIONS OF MOTION 

The six degree of freedom equations of motion for the submersible vehicle shown in 
Appendix A were taken from Smith, Crane, and Summey [Reference 1:11-16]. 
Differentiation with respect to time is denoted by a dot over the variable; e.g. u = 

These equations are referenced to a right-hand orthogonal axis system fixed in the body 
(vehicle) as shown in Figure 1 . Since these equations are in reference to a body fixed axis 
system, the Euler angles of pitch (0),roll (<[)), and yaw (^) are used to specify orientation 
with respect to the inertial reference system. The rotation sequence for 9 , 6 and and the 
Euler angle rates for 9, 0 and shown in Appendix B were taken from Smith, Crane, and 
Summey [Reference 1:18-20]. Major variables and parameters as defined by Smith, 
Crane, and Summey [Reference 1:7-10] arc given below : 

1. Dynamic Variables 

u,v,w - Linear velocity components of vehicle with respect to origin of 

body a.xes system relative to fluid. 

p,q,r - Angular velocity components of vehicle with respect to body 

axes system relative to inertial reference system. 

X,Y,Z - Hydrodynamic force components along body axes. 

K,M,N - Hydrodynamic moment components along body axes. 
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2. Mass Distribution Parameters 



m 


- Mass of ihe Hooded vehicle, including the mass of the 
entrained fluid. 


W 


- Weight of the flooded vehicle, including the weight of the 
entrained fluid ( = g m ; where g Is the acceleration of gravity). 


V 


- Displacement volume of the vehicle. 


B 


- Buoyancy force acting on the vehicle ( = p g V ). This is 
independent of the inertial mass distribution of the submersible 
vehicle, including whether or not it is flooded. 




- Coordinates of the CG (center of gravity) in the body axis 

system (Figure 1). These will depend on the mass 
distribution of the vehicle, including the mass of the entrained 
fluid. 




- Coordinates of the CB (center of buoyancy) in the body axis 

system (Figure 1). These are independent of the mass 
distribution system, but may vary with the addition or removal 
of external appendages. 


I ,I ,I 

x’ y’ z 


- Moments of inertia about the body system axes, including the 
entrained fluid. 


I ,i ,i 

xy X/ yz 


- Products of inertia about the body system axes, including the 
entrained fluid. 



3. Remaining Parameters 



P 


- Mass density of fluid medium 


1 


- Reference length used to nondimensionalize the hydrodynamic 
coefficients. 



b(x), h(x) 



^nose’^tail 



- Width and height of vehicle in its xy and xz planes, 
respectively, at location x measured in the body axes system 
(Figure 1). These quantities are required in the integrals 
defining the crossflow forces and moments in the equations of 
motion, and are tabulated within the Steady State Computer 
Program (Appendix C). 

- Sternplane, bowplane and rudder deflection angles in radians 
(Figure 1). 
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Figure 1. Scheinutic Indicating Positive Directions of Axes, Angles, V'elocities, Forces, and Moments 



II. SYSTEM SOLUTIONS IN THE VERTICAL PLANE 



A. GENERAL 

In the steady state condition, the submersible will have reached constant linear and 
angular velocities. Therefore, the body fixed linear accelerations (li , v, w) and the body 
fixed angular accelerations (p , q , f) will be zero. Similarly, the vehicle will have reached 
a constant angle of pitch (0) making its derivatives (0) equal to zero. Since this analysis is 
restricted to steady state solutions in the vertical plane, the angle of roll (<j)) and its 
derivative (<[>) will be zero (in Chapter IV, the case where the angle of roll (<{)) is 180 
degrees will be discussed during the numerical integration analysis). The angle of yaw (^) 
and its derivative (^’) will likewise be zero due to the vertical plane restriction. It should be 
noted that had this analysis not been restricted to the vertical plane, steady state yaw (^) 
would not necessarily be zero, thereby allowing the angular velocities (p, q, r) to be non- 
zero. However, since this analysis was restricted to the cases where <{), 0 , , p, q and r 

are all zero, the equations of motion for six degrees of freedom for the steady state 
condition can be reduced to: 

* Longitudinal (Surge) Equation of Motion: 

(W . B) sin S = ; Xvv v2 + w2 + X^j^ uv6, + uw ( X^j^ 6^ + 6b)] 

* y ( 5^6 s6s * Xab6b «b^ + + ^^Xp,„p 

• Lateral (Sway) Equation of Motion: 

- (W-B) cos 0 sin<t>= Yy uv + Yy^y vw + Y^j. u2 6[-j 
f^nose 

- 1 [ CDy h(x) (v)2 +Cj)z b(x) (w)2 , dx 

•'xtail 
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Normal (Heave) Equation of Motion 

- (W-B) cos 0 cos ({) = uw + Zyy v2 + u2 (z^^, 
pnose ^ 

L ^Dy h(x) (v)2 +Cd 2 b(x) (w)2 - dx 

ixiaii 

Roll Equation of Motion: 

- (yG W - yg B) cos 0 cos (j) + (zQ^ ' ^B^) cos 0 sin ^ 
= ^Kv uv + Kvw vwj + u 2 Kprop 

Pitch Equation of Motion: 

(xQ W - xg B) cos 0 cos (j) + (zQ W - zg B) sin 0 = 

[ Mw uw + Mvv v 2 + u 2 (m^^ 65 + M^j^ 65 ): 

J '*^nose 

Cgy h(x) (v)2 +Cgz b(x) (w)2 ] x dx 

^tail 



Yaw Equation of Motion: 

- (xQ W - xg B) cos 0 sin <{) + (yQ W - yg B) sin 0 = 
Ny uv + Nyw vw + N^j. u2 6r + u2 Nprop 



mosc , 



^^tail 



Cgy h(x) (v)2 + Cdz b(x) (w)2 I — x dx 



B. CONDITIONS 

1. Defining Additional Terms 
a. Excess Buoyancy, bB 

Excess buoyancy is defined as 6B = B - W where B is the submersible's 
total buoyancy and W is the submersible's total weight. 
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b. Longitudunal Center of Buoyancy, xq^ 

The longitudinal center of buoyancy is defined as XQg ~ '''here 

xq is the longitudinal center of gravity with respect to the body fixed axis and Xg is the 

longitudinal center of buoyancy with respect to the body fixed axis. 

c. Vertical Center of Buoyancy, 

The vertical center of buoyancy is defined as ZQg where Zq is 

the longitudinal center of gravity with respect to the body fixed axis and Zg is the 
longitudinal center of buoyancy with respect to the body fixed axis (zQg is assumed to be 

positive). 

2 . Assumed Conditions 

a. Lateral Centers of Gravity, y^ , and Buoyancy, y^ 

The lateral center of gravity and the center of buoyancy are assumed to be 
on the same centerline plane (y^ = Yg = 0). 

b. Propeller Speed,n (revolutions per minute) 

The propeller speed is assumed to be zero (n = 0). 

c. Propeller Coefficients, Kp^^p and Np^^p 

From Smith, Crane, and Summey [Reference 1:30], the propeller 
coefficients are zero (Kp^op = Np^op = 0). 

C. REVISING THE EQUATIONS OF MOTION 

Using the term for vertical center of buoyancy (zqb^’ expression (zqW - ZgB) 

may be written as (zQgW - zg6B). Similarly, using the term for longitudinal center of 
buoyancy (xQg) the expression (xqW - XgB) may be written as (xQgW - Xg6B). Also, 
the term u^Xpj-Qp may be written as: 

u^Xprop = u2Coo(Tl2-l)=u2CDo[(^^^^|^^^f-l] = CD0A2n2 - Cdou2 



8 



where A is a constant. 



Since the shaft speed (n) is zero, the expression may be further reduced to 
*J^^prop = ■ Cdou-. Substituting these expressions plus the term for excess buoyancy 

(6B) and the assumed conditions revises the equations of motion for the six degree of 
freedom system as follows: 



Longitudinal (Surge) Equation of Motion: 

-6B sin 0 = ! Xw v 2 + X^w + X^^^ uvb^ + uw ( 65 + X^^j^ 6^) 



4 * 



(Xfc6sV 



+ ^("ibdb *b‘ + \,tit '1: ' Coo u- 



Lateral (Sway) Equation of Motion: 

6B cos 0 sin<|)= Yy uv + Yy^ vw + Y^^. u2 b^' 
f^nose 

- 1 L ^Dy ' u V) 



Normal (Heave) Equation of Motion: 

6B cos 0 cos ([) = Zy,’ uw + Zyyv2+u2(z^^ ^s+Z5h 6b) 



J '-’^nose 

xtail 



^ Cpy h(x) v 2 +Cd 2 b(x) w 2 ! - -W dx 
^ Uci(x) 



Roll Equation of Motion: 

(zQB^ - zb6B) cos 0 sin (p= Ky uv + Ky^y vw] 

Pitch Equation of Motion; 

(xqb ^ ■ ^B cos 0 cos (j) + (zQ3 W - zg 6B) sin 0 = 

[ Mw uw + Myy v 2 + u 2 5 s + M^,^ 6h) 

J ^^nosc ^ 

^ CDy h(x) v2 +Cdz b(x) w2 x dx 

^tail 
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• Yaw Equation of Motion: 

(-XGB W + xg 6B) eos 0 sin <[) = |^ Ny uv + Ny^y vw + u2 6,-] 

J ^^nose . 

■ Coy h(x) v2 + Cd 2 b(x) w2 ! — X dx 

xtail 

These six equations only have five unknowns: u,v,w,0, and (p. Therefore, two of 
these equations must be dependent and additional eonditions are required in order to make 
the number of equations equal the number of unknowns. 

D. ADDITIONAL CONDITIONS 

The next eondition to be applied to the vehicle is that the rudder will remain 
eenterlined, that is 6r = 0. Sinee the solutions of interest are those in whieh the vehicle 
remains within the vertical plane, it can be further speeified that the linear velocity in the 
transverse direction (v) equals zero. Recalling that the angle of roll (<[)) has been previously 
assumed to equal zero, the trigonometrie funetions of <|) ean be identified as sin<t) equals zero 
and eos<{> equals unity. Substituting these quantities back into the equations of motion, 

three of the six equations of motion (sway, roll, and yaw) yield trivial solutions. In 
addition, the eross-flow veloeity term (U^f) for the heave and pitch equations can be 

red need to: 

Ucf(x) = ^(v +xr)- + (w - xq)- = |w| 
since v, r, and q are zero. Furthermore, since is constant, it can be taken outside the 

integral. Therefore, the three remaining equations ean be wriiten as: 

• Longitudinal (Surge) Equation of Motion: 

-6Bsi„eJ '‘7 ’'was 

. + ( ^6s6s + ^6bbb - ^do 



10 



Norniiil (Heave) Equation of Motion: 

6B cos B = uw + u 2 65) - Cp 



'^nose 



z I II Jx 

w 

^tail 



Pitch Equation of Motion: 

(xGB W - XB 6B) COS 0 + (zGB W - ZB 6 b) sin 0 = 

/ ^nosc 

Mw uw + u2 (m^^ 6s + M^j^ 6b)+ Cq/ I b(x) x dx 

/X.oil 



E. VERTICAL PLANE EQUATIONS OF MOTION 






By defining the terms A^ as the | b(x) dx and x^ as b(x) x dx, the three 

' Xfail 

remaining equations defining motion in the vertical plane can be written as: 



• Longitudinal (Surge) Equation of Motion: 

- 6B sin e = Xww w 2 + uw ( 65 + 65) 

+ ““ ( ^6s6s + >^6b6b ^b') - Cdo U“ 



Normal (Heave) Equation of Motion: 

6B cos 0 = Zw uw + u 2 (z^^ b^+Z^j^ 6^) - w |wl A,v 



• Pitch Equation of Motion: 

(xgb - XB bB) COS 0 + (a3bW-zb6b) sin 0 = 

Mw uw + u2 (m^^ 6s + M^J^ 65)+ Cd^ w |w| xaA,,. 

F. COMPUTER PROGRAM DEVELOPMENT 

Taking these three equations which describe motion in the vertical plane, solving the 
first two for sine 0 and cosine 0, respectively; and dividing all three through by u- yields: 
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Longitudinal (Surge) Equation of Motion: 

e _ j_ 1 6s + 6b) ■ 

u2 6B [ + ( 6s^ + 6b^) - Cixi 



Normal (Heave) Equation of Motion: 



cos 6 _ 1 
u2 6B 






w 



(Z6s 



'6b 



^b) -Cd2 



W W 



Av 



• Pitch Equation of Motion: 

(xQB W - xg 6 b) + (zQB W ■ zg ^ = 

u^ u^ 

M,v + (m^^^ 6s + M^b ^b)+ Cdz xaK- 

Now, defining the quantity ^ as w', and substituting w' back into the three 



equations: 



if w is positive; 



• Longitudinal (Surge) Equation of Motion: 

sin e = . JL I ^w6s + ^w6b ^b) 

u2 6B L +(^6s6s^s^+ ^6b6b ^b^) ■ ^DO _ 

• Normal (Heave) Equation of Motion: 

= L'Av(w') + (Zj,j 5s+Z(,^ 5h) -Cdz(w')2aJ 
U“ 6B 

• Pitch Equation of Motion: 

(xqb W - xg 6 b) £2S-S- + (zQg W - zg 6 b) -^h-S- = 

u2 u2 

Mw (w’) + (m^^ 6s + M^J, 6b) + Cg 2 (w')2 xaAw 

however, if w is negative: 



Longitudinal (Surge) Equation of Motion: 
sin8_X ' Xww(-')2+(W)(x^j^5s + X„5^6i,) I 
6B ' +(x5s5s6s^+ Xj|,j(, 6b^) ■ Cdo j 
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Normal (Heave) Equation of Motion: 

“S0 = 1 Zw(w') + (Z^ <>s+Z(,^ 6|,) +CDz(wr A., 



Pitch Equation of Motion: 

(xGB W - xg 5 b) + (^GB " ^B = 









M,v(w')+ ftb)-CDz (w')2xaA,v 



yn 0 cos 0 . 

Substituting the equations for ^rid ^ 2 '^ito the pitch equation yields the 



u- 



following expressions: 



if v> is positive: 

(xgb W - XB dB) I Zw (w') + (Z^sjj. ds+Z^j^ d^) - Cdz(w')“ A,,, 



dB 



. (zgbw.zb.b) -L 



dB 



( ^dsds 



+ ^dbdb^b 



Mw (w') + (m^j, ds + dj + Cdz (W)^ xaA,v 



and if w is negative: 

(xGBW-XBdB) 1 Z\v (W) + (Z^s ds+Z^j^ db) + Cdz(W )-A« 



dB 



w 1 ^ww(w')2+(w')(x^^^^ds + X^^,^j,db) 

IzGB W - ZB dB) ^ ^ 2 x 2\ ^ 

dB + ( XacAc ds + A^j^bb "b ) ■ ^IX) 



■"l^dsds'N " 
Mw(w')+ (M^, 6i + M(,j, i>„)-CDz (wf XaA«. 



Rearranging these two equations to get them into the form: A(\v')- + B(w') + C = 0. the 



expressions become: 
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if w is positive: 

xaAv, + (zgb W - ZB 6B) -L X 

6B 

+ (xgb W - XB 6B) -1 - Cdz A.V 
6B 



ww 



(w'f 



r Mw + (zgb W - zb 6B) ^ 65 + 6b) 



6B 



(xgb W - XB 6b) J^Zvv 



(w') 



L 6B 

(zgb W - ZB 6B) X53J,, 5,2 + 61,2 . C„„) 



6B 



+ (Mbs 6b) - (xgb W - XB 6B) ^ (z^^, 6b) 

6B j 



and if w is negative: 

I - Cdz xaAw + (zgb W - zb ^b) -L- Xww 

' 6B (w')^ 

I - (xgb W - XB 6 b) Cdz Av 

L 6B J 

Mw + (zgb W - zb 6b) ^ ( X^5^ 6s + X^^^ ^b) ] 

+ 6B |(w'i) 

- (xgb W - XB 6B) j j 

6B 



= 0 



(zgb W ■ ZB 8B) M Xj^j, 6,2 + 6i,2 ■ Cu„) 

oB I 

+ (Mj, 6s+Mj,^ 6[,).(xgbW.xb 6B) 6,+Z;,|, 6b')’ 



= 0 



6B 

These quadratic expressions were then solved using the equation: 
. -B ± \ B 2 - 4 AC 



w = 



2 A 



where: 



B =6B M„ + (zGB W - zg 6B) ( X^j^ 6, + ^b) 

-(xgb W -xb6b) Z* 

C = (zgb W - ZB 6B) ( Xf,j.5, 5,2 + X^bb,, 6b^ - Coo) 

+ 5 B (Mj^ 5 , Mjb 6b) - (zqB W - XB 5 B) (z^, Sj+Z^b 6b) 

if w is positive: 

A = 6B(Cdz xaAw) + (zgb W - zb 6b) X^w + (xGB W - xb 6b) C^z Av; 
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and if w is negative: 

A = - 6B(Cd^ xaA,v) + (zGB W - 5 B) X^av * (xGB ^ ' ^B ^Dz 

The value ol'B was determined using the computed values of w' and the equation for 



tangent 6: 



sin 0 

tan 0 = 

cos 0 



II I . cos 0 

How'ever, the value ot 

u- 



varied depending on the value of w', which lead to two 



possible solutions: 

Equation for tan 0 if w’ is Positive: 



Xw„(w'|- 


-(w' 


l( X„6s 


^S + ^'wbb ^b) 


'■( 


.^bsbsi'j 


'> 

s" + 


^bbbb 


bb^) + Coo 




Z^v 


(w')+ + 




^b! 


1 -CDz 


Aw(w')|w'| 




Equation 


for tan 0 if 


w' is Negati 


ve: 










Xww(W)2 


- (w' 


'!.Xvvft,s 


^ wbb ^b) 


i-J 


^bsbs ^ 




^bbbb 


bb~) + Coo 




Z\v 


(w') + + 


(^bs ^s+Z^,^ 


bb) 


+ Cdz 


Av(\v'|w'| 





In either case, the value of u^ was computed using the expression derived from the 
surge equation of motion; 



u- = 



5B sin 0 



- Xww (w')2 - (w')( X^,^^ 65) - ( bj- + 6i,“) + Coo 



This leads to tw'o possible solutions for u (i.e. u = ± \u-). The value of w was 
computed using w = u (w'). Q)mbining the two possible solutions of u with the two 
possible values for w' derived from the quadratic expression lead to four possible 
combinations of solutions for u and w. The computer program which calculates these four 
possible solutions is contained in Appendix C. It is an interactive program designed to 
allow the operator to select the amount of excess buoyancy as a percentage of vehicle 
weight (6B), the deOection of dive planes in degrees ((\), the ratio of bow planes to dive 
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weight ( 6 B), the deflection of dive planes in degrees (h^), the ratio of bow planes to dive 
planes ( 65 / 65 ,), the location of x^g and Xg from body fixed axis origin as a percentage of 
length, and the location of z^g and Zg from body fixed axis origin in feet. 

G. STEADY STATE RESULTS 

Figures 2, 3, and 4 show typical steady state solutions for surge velocity, heave 
velocity, and pitch angle as a function of dive plane angle. The two cases shown vary the 
location of the longitudinal center of buoyancy; for case (a): x^^g = - 1 % of the vehicle 

length (L), and for case (b): x^g = + 1 % L. The following parameters were the same for 

both cases; excess buoyancy , 6 B = 2 % of the vehicle weight (W); deflection of bow 
planes, 65 = 0 ; location of horizontal and vertical centers of buoyancy. Xg = Zg = 0 ; and 
location of vertical center of gravity, z^^g = 0.1 feet. 

All runs developed four solutions. For two of the solutions, the magnitude of the 
surge vclocitities were large while the magnitudes of the associated heave velocities were 
relatively small. This has been descibed by Booth [Ref 2: 297J as "predominantly forward 
motion". The other two solutions had small surge velocity magnitudes and large heave 
velocity magnitudes. Booth [Ref 3: 346] referred to this type of motion as "nearly vertical 
ascents". The positive or negative nature of the velocities are associated with the value of 
pitch angle. Positive heave velocities are associated with pitch angles greater than 90 
degrees. That is the submersible would be ascending in a belly up orientation. Although 
this steady state analysis computes four possible solutions, it gives no indication as to 
which of the solutions are stable (if any). Dynamic response and stability criteria will be 
discussed in the next chapter. 
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Figure 2. 



Steady State Vertical Plane Solutions for Surge Velocity (u) 
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Figure 3. Steady State Vertical Plane Solutions for Heave Velocity (w) 
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Figure 4. Steady State Vertical Plane Solutions for Pitch Angle (0) 



19 



III. DYNAMIC STABILITY 



A. GENERAL 

The first portion of this chapter uses the six degree of freedom equations of motion 
along with the Euler angle rate equations for the derivatives of the angles of pitch and roll 
(0 and <[)) to predict the dynamic stability when movement is restricted to the vertical plane. 
These equations of motion are then linearized around the vertical plane steady state nominal 
points computed in chapter II. Eigen value analysis is then used to compute the stability of 
each solution. The second part of the chapter expands this analysis to include all six 
degrees of freedom and uses the same steady state nominal points to predict the dynamic 
stability when motion is not restricted to the vertical plane. 

B. RESTRICTING MOTION TO THE VERTICAL PLANE 

Since motion is restricted to the vertical plane, the body-fixed transverse veloctiy (v) 
and its derivative (v) are zero. The angles of roll (<{)) and yaw (V ), and their derivatives (<j) 
and are also zero. We will continue to assume that the lateral center of gravity and the 
lateral center of buoyancy arc on the same centerline plane (y^ = yg= 0), and the rudder is 
centerlined (hr =0). 

C. LINEARIZED VERTICAL PLANE EQUATIONS OF MOTION 
Substituting these values into the original equations of motion (Appendix A), yields 

trivial results for three of the six equations: Lateral (sway), Roll, and Yaw. The remaining 
equations reduce to the following form: 
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Lxingitudinal (Surge) Equation of Motion: 

[ m - Xu u + m zq q = [ -Coo + +X^o6b^l^^ ^2 

+ + X,,, 5 o‘'>h uw + X^^,6s + uq+ [X^^v. w2 

+ Xwq - m wq-i- Xqq + mxo q2 - (W - B) sin 0 

Normal (Heave) Equation of Motion: 

r 

m - Zvv w - m xg + Zq q = u^ 

+ Z^v uw + m uq+ mzQ q- 

^xnose 

- 1 lQ)z b(x)(w-xq)- dx + (W - B) cos 0 

✓ Mail 



Piteh Equation of Motion: 

mzo u - Mw + m xg w + ; ly - Mq q = ; M^v^^bs + M^^db u^-i- M^ 



+ Mq - mxG uq- mzG 




Coz b(x)(vv-xq)2 



(w ^xq) 

Ud<x) 



X dx 



uw 



- (xgW - xbB) eos 0 - (z{,W - zbB) sin 0 



These three equations which describe motion restricted to the vertical plane are functions of 
four variables (u,w,q,0) plus their deriatives (u,w,q,0). The equations of motion and the 
Euler angle rate equation for 0 were linearized using the following generalized proeedure: 



Bll(i,l)u + Bll(i,2)w + Bll(i,3)q + B1 l(i,4)0 




i f'fi ^ 
u + !^ 

Uo JV . U 



W + 



^fi 

0q 



q + 

q» 



e 

00 



where the B1 l(i)'s are the constants associated with the derivatives of the variables, the 
functions fj represents the right hand side of the nonlinear equations, and i = 1 to 4 
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identifies the three equations of motion (surge, heave, and yaw) plus the Euler angle rate 
equation for 0. The partial derivatives were computed as follows: 



Partial Derivatives of the Lxtngitiidinal (Surge) EOM: 

- ^ = 2(-Coo + ^6s6s^^ +^6b6b^t* W + (^w6s^^ + ^wbb^^Wo = A1 1(1,1) 
du 

= 2X,^wo + (X„^,6s + X,.5b6b)uo = A11(J,2) 

dw 

= (>^wq - m)wo+ (Xq^^s^s + Xq^jb^bjuo = A1 1(1,3) 

— ^ = (W - B) cos 00 = A1 1(1,4) 

60 

Partial Derivatives of the Normal (Heave) EOM: 

= Z,.vvo + 2(z,,6s + Z,b6b)uo = All(2,l) 
du 

1^2 = uo4v - 2 CdAv|wo| =A1 1(2,2) 
dw 

= uo(Zq + m) + 2 CdzAwXa1wo| = A1 1(2,3) 
dq 

=-(W-B)sin 00 = A1 1(2,4) 

60 

Note on differentiation procedure: The cross-flow velocity term (U^f) was 
reduced to Uei(x) = ■ (v +xr)- + (w - xqp/*'^ = [w^j®-^ = |w - q| . Allowing 
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the integral term in the heave equation to be defined as 1,, where; 



Ii=Cnzj b(x)(w-xq)|w - 

J xtail 


=2CdzI 

d(w - q) J 


^xnosc 

! b(x)o 

Mail 


Ah = 


1 1 


j d(w_- q)l 


dw 


\ 6(w - q) 1 


1 dw 1 


Ah- ^ 




i 6(w - q)\ 


dq 


d(w-q) /I 


1 dq ) 



f 

¥ xt-j 



J ^xnose 

b(x) dx = 2 Cdz |\V()| A„ 

Mail 



* Partial Derivatives of the Pitch EOM: 

= 2 (m^,5s + M,i,6b)uo + M, W(, = A1 1(3,1) 

= M,v uo +2 uoCdzA.vXa!woI = A1 1(3.2) 
dw 

= (M^ - mxo)u(, -(m/G)wo - 2Ciy, 1a!wo| = A1 1(3,3) 

= (xqW - xnB) sin0o - (zc.W - zbB) cos 0o = A1 1(3,4) 

60 

Note on differentiation procedure: Again the cross-llow velocity term (U^f) was 
reduced to |w - q|. Allowing the integral term in the pitch equation to be defined 
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as where: 

J M xnosc 

b(x)(w-xq)lw - xql x dx 

xt: 



ai 







/ 


<^W 


\ a(w - q) 


dh ^ 


\ 


aq 


a(w - q) / 



I 



‘I 



= 2CdJwo 1 b(x) X dx = 2 Cdz Kol xa Aw 



‘1 



= - 2CdjJwo1 x 2 b(x) dx = - 2 Cdz kol U 



• Partial Derivatives of the Euler Angle Rate Equation for 0: 

= 0 =A1 1(4,1) 



^^4=0= Al 1(4,2) 

aw ^ ^ 



a £4 
aq 



1 = Al 1(4,3) 



=0 =A1 1(4,4) 

ae 



The constants corresponding to the derivatives of the four variables (i.e. the left hand 
side of the four equations) are as follows: 
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• Constants Associated with Derivatives for the l>ongitiidinal 
(Surge) EOM: 



li ^ m -Xu = Bll(l,l) 



q => mzG = Bl 1(1,3) 

B1 1(1,2) = B1 1(1,4) = 0 

• Constants Associated with Derivatives for the Normal 
(Heave) EOM: 



w => m - Z^. = Bl 1(2,2) 
q => - (mxG + Zq) = B1 1(2,3) 



B1 1(2,1) = B1 1(2,4) = 0 



• Constants Associated with Derivatives for the Pitch EOM: 



li => mzG = B1 1(3,1) 



w => - (M\v + m xg) = B 1 1 (3,2) 
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ci => Iy-H^ = B1 1(3.3) 
B 11(3,4) = 0 



• Constants Associated with Derivatives for the Euler Angle Rate 
Equation for 6 : 

0 => 1 = B1 1(4,4) 

B1 1(4,1) = B1 1(4,2) = B1 1(4,3) =0 

These expressions can be arranged in a matrix format to form the linearized equations 
of motion in the vertical plane about the nominal steady state points. The matrix format is 
as follows: 

Bll X X 1 = All X XI 

where: 



All(l,l) 


A1 1(1,2) 


A1 1(1,3) 


All(l,4) 


A1 1(2,1) 


A1 1(2,2) 


A1 1(2,3) 


A1 1(2,4) 


A1 1(3,1) 


A 11(3,2) 


A1 1(3,3) 


A1 1(3,4) 


A1 1(4,1) 


A1 1(4,2) 


A1 1(4,3) 


A1 1(4,4) 


Bll(U) 


Bll(l,2) 


B1 1(1,3) 


B1 1(1,4) 


Bll(2,l) 


B1 1(2,2) 


Bll(2,.3) 


B 11(2,4) 


B1 1(3,1) 


B1 1(3,2) 


B 11(3,3) 


B1 1(3,4) 


Bll(4,l) 


B 11(4,2) 


Bll(4,.3) 


Bll (4,4) 
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and 

u 

X]= 

‘I 

0 

D. VERTICAL PLANE COMPUTER PROGRAM DEVELOPMENT 

The matrix format of the linearized equations of motion was used in the computer 
program shown in Appendix D to predict the dynamic stability of the vehicle with motion 
restricted to the vertical plane. The program is interactive in that it allows the operator to 
select which of the four data files from the steady state analysis (Chapter II) will be used to 
define the nominal points for the linearization process. An eigen system subroutine was 
used to find eigen values and eigen vectors. The program's output was called the degree of 
stability and only shows the largest real part of all eigen values. The stability criteria is 
sueh that the degree of stability must be negative in order for the solution to be stable. 

E. VERTICAL PLANE DYNAMIC RESULTS 

Of the four possible steady state solutions computed in Chapter II. only one solution 
from each case yielded stable characteristics. There were some cases in which none of the 
solutions were stable for certain ranges of parameters. The general trend of the linearized 
dynamic results are fairly well demonstrated by the two cases diseussed in Chapter II. 
Recalling the parameters of these cases: 6 B = 2%, ratio of bow to stern planes = O. 

Zgb = 0-1 f^ct, and ^3 = ^3 = The first ease placed the longitudinal center of gravity aft 
of the longitudinal eenter of buoyaney (Xqb = - 1 %), and the seeond ease plaeed the 
longitudinal center of gravity forward (x^b = + 1 %). Once again dive plane deflection 
angle ( 6 ^) was varied from - 20 to + 20 degrees for both cases. Figure 5 shows 
longitudinal velocity (u) as a function of dive plane angle (5^). Case One (Xqb = - I %) 
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showed predominantly forward motion, while case two (x^g = - 1 %) yielded nearly 

vertical ascents. Figure 6 shows vertical motion (w) as as a function of dive plane angle 
(6^). The results concur with Figure 5, case one shows very little vertical motion while 

case two demonstrates a larger value. It is interesting to note that vertical motion (w) for 
case one is positive for dive plane angles between - 20 and -4 degrees. The case one values 
of 9 shown in figure 7 for dive plane angles between - 20 and -4 degrees concur with this 
observation. The values of 0 greater than 90 degrees indicate the submersible is ascending 

in the belly up position. The stability in the vertical plane is shown in figure 8. Degree of 
stability (e) is shown as a function of dive plane angle (6^). This figure shows both cases 

to be stable in the vertical plane. 




Figure 5. Stable Vertical Plane Solutions for Surge Velocity 
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Figure 6. Stable Vertical Plane Solutions for Heave Velocity 
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Figure 7. Stable Vertical Plane Solutions for Pitch Angle (6) 
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Figure 8. Degree of Stability in the Vertical Plane 

F. LINEARIZ.ATION OF FULL EQUATIONS OF MOTION 

Referring back to the complete equations of motions shown in Appendix A, these six 
equations which describe motion for the six degree of freedom system are functions of 
eight variables (u,v.w,p.q.r,<I),0) plus their deriatives (u,v,w,p,q,f,<}>,0). The six 
equations of motion (Appendix A) along with the Euler angle rates (Appendix B) were 
linearized as folows: 

bjlii + bj2v + bj3w + bj4p + bj5q + bj6f + bj7<}) + bj80 = 



where the bj's are the constants associated with the derivatives of the variables, the 
functions gj represents the right hand side of the nonlinear equations, and j = 1 to 8 
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identifies the six equations of motion plus the Euler angle rate equations for <|) and 0. 
respectively. The partial derivatives of these equations were computed as follows: 



Partial Derivatives of Longitudinal (Surge) EOM: 



=(^w6s + ^\v6b f'hVo ds + 5b )u() - ZCdouo = 

ou 



= 2X,v•^vWo + 5s + X^,^h = al3 

(7 W 



= -mwo + X.vqWo + (x >,5 5s + X^^ 5b)uo = al5 



dq 



^gi 

d 0 



= - (W - B) cos 0 = al8 



dv dp 



= 0= al4 



dgl 

dr 



= 0 = al6 



dgi 

d(p 






• Partial Derivatives of Lateral (Sway) EOM: 



=YvUo + YvwW'(, - Cd.Aw |w„| = a22 
dv 



=mwo + Yp uo +Y,vpW(, = a24 

dp ' ' 



= - muo +Yr U()+ Yvvr wq - Cl)^AwX,^!wo!= a26 
dr 



= (W - B) COS0 = a27 
d(J) 



= al7 
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^g2 

du 



= 0 = a21 



M-n- 



dw 



= 0 = a23 



= 0 = a25 
dq 



%2 

d0 



= 0 = a28 



Note on dilTcrcntiation procedure: The cross-flow velocity term (U^j.) is given 
by: Ucf(x) = [(v +xr)- -t- (w - The integral term in the sway equation 



was defined as I^, where: 



1.2 = 



Cdv h(x) (v+xr)2 +Cd^ b(x) (w-xq)2 ; ^ dx 



J^^nose 

xtail 

J '^nose 

xtaii 



- Ud(x) 



dl^ 


( -'M 


/ V + xr] 


1+ I ^ ( 


' V + xrA 


dv 


\dv / 


i Ucf i 


' dv 1 


1 Ucf / 



= Cdz Av\ Wo“ 



= CdzW() 2 Aw = Cdz |wo| A„ 

W(,-^ 



Ucf - (v - 1 - xr) — Ucf 
dv i 

Ucf ^ 



51.2 _ 


1 


( V + _xr j 




' V + xr\ 


dr 


\ dr j 


i Ucf ] 


dr ' 


i Ucf / 



= Cdz AwWo^ ^ 



u 



cf 



.dJJcf^ 

’ dr 



Wo^ 



dl, 



dw \ dw 



di 



' A( ^ + xr) \ _^^ 



Ucf / dw \ Ucf 



= f 1 f ^ + xQ ] ^ 0 

dq 1 dq / 1 Ucf ) dq \ Ucf 



Partial Derivatives of Normal (Heave) EOM: 



= Zvv Wo + 2 (Zos 6s + 6b)u(, = a31 

du 

= 4vUo - 2CDzA,v|w'o| = a33 



= m uo + Zo U() + 2 Coz A,v xa!w(,! = a35 
dq ^ 

(w - B) sin 00 = a38 
(36 



^ Q ^ 332 f'g3 = 0 = a34 
f)v c)p 



''S3 ^ 0 = a36 
dr 



= 0 



Note on dilTercntialion procedure: The integral term in the heave equation 
was defined as where: 



U = 



[ 



^nosc 



Cdv h(x) (v+xr)2 +Cdz b(x) (w-xq)2 , dx 



xtail 

J '^nose 

xtail 



■ Ucltx) 






dl4 

dv 



( il W(w.xq« , X- 

\ dv / \ Ucf / dv 1 Ucf / dv 

= - Coz AvvW'o^ — Wo = 0 because — = 0 



Wo*- 



dl4 ^ 


=( d) 




f \v - xq \ 




\ aw / 


\ Ucf 1 dw \ 


i Ucf ) 



= 2 Cdz A^vv - xq)'^ + I ^ I Ucf - (w - xqf 

Ucf Ucr L 

= 2 C[)^ A^.wo2-J- +Cdz A^.vvo2 -J J |vv(,| - wq'^^ ^ = 2 Cdz A^|wo| 

Fo! wo2l Fo! 

because = ^Ucf{x) = (v +xr)^ + (\v - xq)^/’'^2(w - xq) = |-'| 

= Cdz AwWq 2 J - X Ucf - (w - xq}^- 

Uef / Ucf 

= Cdz AwW()2^'aJ^^ = Cdz AcvXa|W()1 



du 



dl 

dq 



w - xq 



+ 1 



dq 



il4 

dr 



= 0 



• Partial Derivatives of Roll EOM: 



= KvUo + KvvvWo = a 42 
dv 

= - m zg W() + Kp uo + Kwp Wo = a44 

dp ^ ' 

d^A 

'r = m zq Uo + Kr uo + K^r wo = a46 



= - ( ZG w - ZB B) COS0O = a47 
d<p 



dg4 

du 



0= a41 



^ = 0 = a43 = 0 = a45 

d\v dq 



= 0 = a48 
d0 
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Partial Derivatives of Pitch EOM: 



= Mw W(, + 2 6s + 6b)u(, = a51 

du 

= Mw uo + 2 Cdz AwXa'woI = a53 
dw 



f-- =- mxG Uo - mzG wo + uq 
dq ^ 



- 2 Cdz IaIwoI = a55 



= (xgW - xbB) sin 0o - (zgW - zbB) cos 0o = a58 
d0 



= 0 = a52 
dv 



= 0 = a54 
dp 



lS5 _ n _ 



dr 



= 0 = a56 



'fe = 0 = aS7 
d(]) 



Note on differentiation procedure: The integral term in the pitch equation 
was defined as where: 



f^nose . _ 

Is = I Coy h(x) (v+xr)2 +Cd^ b(x) (w-xqp ; ^ x dx 

•Ixtail 

J ^^nose 

> 



xtail •'>^tail 



nose 



L X dx 



dl.s _ dl 4 
dw dw 



Xa= 2CdzAwX'a|wo| 



dls _ dl4 
dq dq 



xa = Cdz AwXa^IwoI = Cdz I/^JvvoI 



dl5 _ dl4 
dv dv 



Xa = 0 



dis dl4 „ 
^ — = . — xa= tJ 
dr dr 



35 



Partial Derivatives of Yaw EOM: 



^6 

dv 



Nv U() + Nvw wo - Cdz Aw xa |wo| = a62 



= m xg Wo + Np uo + Nwp wo = a64 

dp F K 



- m XG Uo + Nr uo + Nwr wo - Cdz IaIwoI = a66 
dr 



(xg w - XB B) COS0O =a67 
d(p 






du 



= 0 = a61 



= 0 = a63 
dw 



^ Q ^ 3^5 

dq 



^ = 0 
d6 



Note on differentiation procedure; The integral term in the yaw equation 
was defined as where: 



f^nosc 

U = I i Coy h(x) (v+xr)2 +Cd 2 b(x) (w-xq)2 ] Y_+xl x dx 

•'xtail 



lose f^nose 

il •'^tail 



^^nose 

'xtail 



X dx 



Xa = Cdz |wo! Aw Xa 

dv dv 



f-- = Xa = Cdz AwXa^ iwoi = Cdz IaIwoI 
dr dr 



dw dw 



d^ ^ dh 
dq dq 



xa = 0 



= a68 
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Partial Derivatives of Euler Angle Rate for <1) 

^g7 



dp 

dr 



= I = a74 



= lan B() = a76 



= 0 = a7 1 




dg7 




dv 


dvv 


= 0 = a75 






dq 




dO 



Partial Derivatives (>f Euler Angle Rate for 0 : 



= 1 = a86 
dq 








CC 

II 

O 

11 

X 


= 0 = a82 


'^88 ^ 0 ^ 


f'gs 


du 


dv 


dw 


dq 


^88 = 0 = ag6 


^'88 = 0 = a87 


*^88 ^ 0 = a88 




dr 


d(f) 


d0 





The consianis corresponding lo ihe dcrivaiives of ihe eight variables (i.e. the left hand 
side of the eight equations) are as follows: 
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• Constants Associated with Derivatives for the Longitudinal 
(Surge) EOM: 

u => (m-Xu)=bll q => (mzG) = bl5 

bl2 = bl3 = bl4 = bl6 = bl7 = bl8 =0 



• Constants Associated with Derivatives for the Lateral (Sway) EOM 

V =;. (m - Y;) = b22 p => (-mzo - Yp) = b24 

r => (mxG - Y;.) = b26 
b21 = b23 =b25 = b27 = b28 = 0 

• Constants Associated with Derivatives for the Normal (Heave) 

EOM: 

w => (m - Z^) = b33 q (' * Zq) = b35 

b31 = b32 = b34 = b36 = b37 = b38 

• Constants Associated with Derivatives for the Roll EOM: 

V => (-mzG - Kv) = b42 f ( - Kf) = b46 

p (Ix-Kp)=b44 

b41 = b43 = b45 = b47 = b48 



• Constants Associated with Derivatives for the Pitch EOM: 
li => (mzG)=b51 w => (- mxG - Mvi) = b53 

q (Iy-M,) = b55 

b52 = b54 = b56 = b57 = b58 



• Constants Associated with Derivatives for the Yaw EOM: 

V => (mxG - Nv) = b62 p => ( - Np) = b64 

( => (I^ - NO = b66 

b61 = b63 = b65 = b67 = b68 = 0 
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Constants Ass(»ciated with Derivatives Idr the Euler Angle Rate 



Equation for <}> : 

(p 1 = b77 

b71 = b72 = b73 = b74 = b75 = b76 = b78 = 0 

• Constants Associated with the Derivatives for the Euler Angle 
Rate Equation for B : 

e 1 = b88 

b81 = b82 =b83 = b83 = b84 = b85 = b86 = b87 = 0 
These expressions ean be arranged in a matrix formal to define the dynamic equations 
of motion for the six degree of freedom system linearized about the nominal steady state 
points. The matrix format is B x X =A x X, where: 



all 


aI2 


al3 


al4 


al5 


al6 


al7 


al8 




a21 


a22 


a23 


a24 


a25 


a26 


a27 


a28 




a31 


a32 


a33 


a34 


a35 


a36 


a37 


CO 




a41 


a42 


a43 


a44 


a 45 


a46 


a47 


a48 




a51 


a52 


a53 


a54 


a55 


a56 


a57 


a58 




a61 


a 62 


a63 


a64 


a 65 


a 66 


a67 


a68 




a71 


a 72 


a 73 


a74 


a 75 


a76 


a77 


a 78 




a8l 


a82 


a83 


a84 


a85 


a86 


a87 


a88 




- 


all 


0 


al3 


0 


al5 


0 


0 


a 18 




0 


a22 


0 


a24 


0 


a26 ; 


a27 


0 




a3I 


0 


a33 


0 


a35 


0 


0 


a: 

cc 




0 


a42 


0 


a44 


0 


a46 i 


a47 


aO 




a51 


0 


a53 


0 


a55 


0 


0 


a58 




0 


a62 


0 


a64 


0 


a66 ; 


a67 


0 




0 


0 


0 


1 


0 


a76 


0 


0 




0 


0 


0 


0 


1 


0 


0 


0 
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bll 


bl2 


bl3 


bl4 


bl5 


bl6 


bl7 


bl8 




b21 


b22 


b23 


b24 


b25 


b26 


b27 


b28 




b31 


b32 


b33 


b34 


b35 


b36 


b37 


b38 




b41 


b42 


b43 


b44 


b45 


b46 


b47 


b48 




b51 


b52 


b53 


b54 


b55 


b56 


b57 


b58 




b61 


b62 


b63 


b64 


b65 


b66 


b67 


b68 




b7I 


b72 


b73 


b74 


b75 


b76 


b77 


b78 




b81 


b82 


b83 


b84 


b85 


b86 


b87 


b88 


- 




bll 


0 


0 


0 


bl5 


0 


0 


0 ^ 

1 




0 


b22 


0 


b24 


0 


b26 


0 


1 

0 ' 




0 


0 


b33 


0 


b35 


0 


0 


0 




0 


b42 


0 


b44 


0 


b46 


0 


0 




b51 


0 


b53 


0 


b55 


0 


0 


0 




0 


b62 


0 


b64 


0 


b66 


0 


0 




0 


0 


0 


0 


0 


0 


1 


0 




0 


0 


0 


0 


0 


0 


0 


1 



xl 


1 

1 


u ^ 




1 

! 1 


V 


x3 




w 1 

i 


x4 




P 


: X5 




q 


x6 




r 


' x7 


1 

1 


4> 


. _x8 ' 


i 


0 
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These matrices may be reordered such that they will be of the form: 



Bll 0 ; ^ XT All 0 X 1 

0 B22 0 A22 X2 



This is accomplished by rewriting the X matrix such that X 1 is the same matix used 
in the vertical plane analysis: 



X = 



xl 

x3 

x5 

x8 

x4 

x7 

x2 

x6 



r u 

w 

q 

e 

■ p 

<t> 

V 



XI 

X2 



The A matrix is restructured into a matrix containing four 4x4 matrices with 
theA12 andA21 matrices containing only the zero clement. 



all 


a 13 


al5 


al8 


0 


0 


0 


0 


a31 


a33 


a35 


a38 


0 


0 


0 


0 


a51 


a53 


a55 


a58 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


a 44 


a47 


a42 


a46 


0 


0 


0 


0 


1 


0 


0 


a76 


0 


0 


0 


0 


a24 


a27 


a 22 


a26 


0 


0 


0 


0 


a64 


a67 


a 62 


a66 
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Similarly, B is restructured into a matrix containing four 4x4 matrices with 



theB12 and B21 matrices containing only the zero clement. 



bll 


0 


bl5 


0 


0 


0 


0 


0 


0 


b33 


b35 


0 


0 


0 


0 


0 


b51 


cr 


b55 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


b44 


0 


b42 


b46 


0 


0 


0 


0 


0 


1 


0 


1 

0 


0 


0 


0 


0 


b24 


0 


b22 


b26 


0 


0 


0 


0 


b64 


0 


b62 


b66 J 



As discussed, the X 1 matrix established to descibe the linearized dynamic stabilty in the 
vertical plane is the same as the XI matrix within X. In addition, the All and Bll 
matrices from the vertical plane analysis are also identical to those in X . That is: 



" All(l,l) 


All (1,2) 


A1 1(1,3) 


A1 1(1,4) ■ 




' a44 


a47 


a42 


a 46 


A1 1(2,1) 


A1 1(2,2) 


A1 1(2,3) 


A1 1(2,4) 




1 


0 


0 


a76 


A1 1(3,1) 


A1 1(3,2) 


All (3,3) 


All (3,4) 




a24 


a27 


a22 


a26 


. All(4,l) 


All (4,2) 


A1 1(4,3) 


A1 1(4,4) _ 




_ a64 


a67 


a62 


a66 _ 



; Bll(l,l) 


B1 1(1,2) 


B1 1(1,3) 


B1 1(1,4) - 




' b44 


0 


b42 


b46 ■ 


B1 1(2,1) 


B1 1(2,2) 


B 11(2,3) 


B1 1(2,4) 




0 


1 


0 


0 


B1 1(3,1) 


B1 1(3,2) 


B1 1(3,3) 


B1 1(3,4) 




b24 


0 


b22 


b26 


B1 1(4,1) 


B1 1(4,2) 


B1 1(4,3) 


B1 1(4,4) 




_ b64 


0 


b62 


b66 J 



The A22, B22, and X2 matrices represent the additional equations necessary to describe 
the linearized motion for all six degrees of freedom; henceforth referred to as the horizontal 
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plane contributions. The eigen I'unetion for the six degree of freedom model is computed 
by taking the determinant as follows: 



det 



All -XBll 
0 



0 

A22-XB22 



= 0 =^(dct All -TBllJIdet A22-TB22 )=0 



Since the eigen function will be the product of these two determinants, the 
resulting eigen values will merely be the union of the vertical plane eigen values and the 
horizontal plane eigen values. The significance of reducing the eigen value calculation from 
an 8 by 8 matrix problem to two 4 by 4 matrix problems is not in the computation time 
saved. But rather in faet that now the horizontal and vertieal stabilities have been separated 
and identified. 



G. COMPUTER PROGRAM DEVELOPMENT 

The matrix format of the linearized dynamie response equations assoeiated with 
motion in the horizontal plane was added to the eomputer program developed previously 
(Appendix D). Onee again, the program is interaetive in that it allows the operator to select 
which of the four data files from the steady state analysis (Chapter II) will be used to 
define the nominal points for the linearization process. An eigen system subroutine was 
used to find eigen values and eigen vectors. Two outputs were added to the program. 
First, the degree of stability in the horizontal plane, and next the degree of stability of both 
planes (i.e., the union of the vertieal and horizontal degrees of stability). Reminder: 
degree of stability must be negative in order for the solution to be stable. 

H. DYNAMIC STABILITY SOLUTIONS 

Q)ntinuing with the same eases from part E, the horizontal stability for the two eases 
arc shown in figure 9. Case two (x^g = + 1 %) is stable in the horizontal plane for all 
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values of dive plane angle ( 65 ). Whereas, ease one (x^g = - 1 %) is unstable in the 
horizontal plane for dive plane angle ( 65 ) between -20 and - 9 degrees. From figure 7, this 
eorresponds to values of 0 greater than 140 degrees. This indieates that the vehiele is 
stable (even in the horizontal plane) for values of 0 greater than 90 degrees. A submersible 
with an angle of piteh greater than 90 degrees will have a negative metaeentrie height and 
will therefore be statically unstable. However, the results shown in figures 7 and 9 indicate 
that the vehicle will remain dynamically stable is such a condition. This 'inverted 
pendulum' type stability will be further investigated during the numerical integration 
analysis (Chapter IV) to see if hydrodynamic and drag forces on the vehicle can actually 
cause this to occur. Figure 10 shows the combined stabilities for the horizontal and vertical 
planes. It should be noted that in general the horizontal plane dictated stabilty for the cases 
considered. 
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degree of slability degree of stnhilUy 



0.15 




delta s 



Figure 9. Degree of Stability in the Horizontal Plane 




Figure 10. Degree of Stability in Both (Horizontal and Vertical) Planes 



45 



Because the real part of the computed eigen values must be negative in order for a 
solution to be stable, it is easy to identify unstable solutions. However, the measure of 
stability for a stable solution is harder to quantify. The magnitudes shown thus far for 
degree of stability seem very small. Does this mean that the solutions are not very stable? 
In an attempt to answer this question, the values for degree of stability for a neutrally 
buoyant submersible (6B = 0 or W = B) were computed. Figure 11 shows the degree of 
stability in (a) the horizontal plane and (b) the vertical plane as a function of surge velocity. 
This figure shows that the magnitude of the values for degree of stability for this neutrally 
buoyant case are indeed of the same order of magnitude as the positively buoyant cases 
discussed earlier. 




0 5 10 15 20 0 5 10 15 20 



Sures Velocity (u) Suree Velocitv (u) 

(a) Horizontal Plane (b) Vertical Plane 

Figure 11. Degree of Stability as a Function of Surge Velocity (u) for a 

Neutrally Buoyant Submersible 
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IV. SIMULATIONS AND DISCUSSION OF RESULTS 



A. DYNAMIC STABILITY ANALYSIS RESULTS 



The dynamic stability analysis included in this Chapter considers stability for both 
planes (i.e. the analysis no longer breaks stability down by horizontal or vertical plane). 
As mentioned in Chapter 111, horizontal plane stability generally dictates the stability of the 
vehicle. 

1. Variations in Longitudinal Center of Gravity 



Figures 12 through 15 show how changing the longitudinal center of gravity 
(x^g) effects the dynamic response of the submersible. For these cases, the amount of 

excess buoyancy (6B) is two percent of wicghl (W), the bow plane deflection angle (6j^) is 
zero, the vertical center of gravity (^i^g) is 0.1 feet, and the longitudinal and vertical centers 
of buoyancy (Xgand Zg) arc zero. The longitudinal center of gravity (x^g) is varied from 
- 2 to + 2 percent of vehicle length. When the longitudinal center of gravity (x^~g) is greater 



than zero, there are stable solutions for the full range of dive plane angles, 
is not true when the longitudinal center of gravity (x^g) is less than zero. 

stable solutions is restricted when x^g= - 0.5 and - 1.0. 



However, this 
The range of 
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(b) \Q^ = 0, - 0.5, - 1.0, - 1.5 and - 2.0 % 

Figure 12. Stable Surge Velocity (u) Solutions for V^ariations in 
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(a) Xq 3 = 0, 0.5, 1.0, 1.5 and 2.0 % 
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delta s 

(b) Xgb = 0, - 0.5, - 1.0, - 1.5 and - 2.0 <7^ 



Figure 13. Stable Heave Velocity (\v) Solutions for Variations in x^g 



49 



theta llieta 



-35 
-40 
-45 
-50 
-55 
-60 
-65 
-70 
-75 

-20 -15 -10 -5 0 5 10 15 20 

delta s 

160 
140 
120 
100 
80 
60 
40 
20 



-20 
-40 

(b) xgb = 0, - 0.5, - 1.0, - 1.5 and • 2.0 % 

Figure 14. Stable .Angle of Pitch (0) Solutions for Variations in 
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Figure 15. Degree of Stablity for Variations in xqb 
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2. Variations in the Amount of Excess Buoyancy (5B) 

Figures 16 through 19 show how changing the amount of excess buoyancy 

(5B) effects the dynamic response of the submersible. For these cases, the bow plane 
deflection angle (5^) is zero, the longitudinal center of gravity (x^g) is - 0.5 percent of 

vehicle length, the vertical center of gravity (zQg) is 0.1 feet, and the longitudinal and 
vertical centers of buoyancy (Xg and Zg) are zero. The amount of excess buoyancy (5B) is 

varied from 0.5 to 2.9 percent of weight (W). The only case where there are stable 
solutions for the full range of dive plane angles is when buoyancy (6B) is 0.5. 
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(a) 6B = 0.5, 1.0, 1.5 and 2.0 
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(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 % 



Figure 16. Stable Surge Velocity (u) Solutions for Variations in 
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(a) 6B = 0.5, 1.0, 1.5 and 2.0 % 
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delta 3 

(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 9r 
Figure 17. Stable Heave Velocity (w) Solutions for Variations in 6B 
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(a) 6B = 0.5, 1.0, 1.5 and 2.0 % 
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(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 % 



Figure 18. Stable .Angle of Pitch (6) Solutions for Variations in 6B 



degree of sij|lulily ^ degree of slnhillfy 




(a) 6B = 0.5, 1.0, 1.5 and 2.0 




(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 % 

Figure 19. Degree of Stablity for Variations in 6B 
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3. \ariations in V'ertical Center of Gravity 

Figures 20 through 23 show how changing the vertical center of gravity (z^g) 

effects the dynamic response of the submersible. For these cases, the amount of excess 
buoyanc} (6B) is two percent of weight (W), the bow plane deflection angle (6^^) is zero, 
the longitudinal center of gravity (x^g) is - 0.5 percent of vehicle length, and the 
longitudinal and vertical centers of buoyancy (Xg and Zg) arc zero. The vertical center of 
gravity (zQg) i.s varied from 0.05 to 0.30 feet. The general trend shown in these figures is 
that the larger \’alues for vertical center of gravity (z^g) have stable solutions over a larger 
range of dive plane angles. 



6 

0.25 
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delta s 

Figure 20. Stable Surge Velocity (u) Solutions for Variations in z^^g 
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Figure 21. Stable Hea\e Velocity (w) Solutions for Variations in 




Figure 22. Stable Angle of Pitch (0) Solutions for Variations in z^g 
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Figure 23. Degree of Stability for Variations in 

4. Variations in the Longitudinal Center of Buoyancy (Xjj) 

Figures 24 through 27 show how changing the longitudinal center of buoyanc)' 
(Xg) effects the dynamic response of the submersible. For these cases, the amount of 
excess buoyancy (bB) is two percent of weight (W). the bow plane deflection angle (6^) is 
zero, the longitudinal center of gravity (x^^g) is - 0.5 percent of vehicle length, the vertical 
center of gravity (z^g) i^' 0.1 feet, and the vertical center of buoyancy (Zg) is zero. The 
longitudinal center of buoyancy (Xg) is varied from - 9 to + 2 percent of vehicle length. 

The general trend shown in these figures is that the positive values for longitudinal center 
of buoyancy (xg) tend to have stable solutions for positive dive plane angles. On the other 

hand, the negative values for longitudinal center of buoyancy (Xg) tend to have stable 
solutions for negative dive plane angles. 
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Figure 24. Stable Surge Velocity (u) Solutions for Variations in Xp 
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Figure 25. Stable Heave Velocity (w) Solutions for Variations in Xg 
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Figure 26. Stable Angle of Pilch (9) Solutions for Variations in 




Figure 27. Degree of Stability for Variations in 
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5. Variations in Bow Planes Deflection Angle (6j^) 

Figures 28 through 31 show how a non-zero bow planes deflection angle (6j^) 

effects the dynamic response of the submersible. For these cases, the amount of excess 
buoyancy (6B) is two percent of wieght (W), the longitudinal center of gravity (x^g) is 
- 0.5 percent of vehicle length, the vertical center of gravity (z^g) is 0.1 feet, and the 
longitudinal and vertical centers of buoyancy (Xg and Zg) are zero. The bow planes are 

given a deflection of - 20 degrees. The significance of the results shown in these figures is 
that for certain dive plane angles (6^ = -3 to -12 degrees) there are two stable solutions. 




Figure 28. Stable Surge Velocity (u) Solutions for a Non-zero Bow Plane 

Angle (6jj = - 20 degrees) 
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Figure 29. Stable Heave Velocity (w) Solutions for a Non-zero Bow Plane 

Angle = - 20 degrees) 
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Figure 30. Stable Pitch Angle (6) Solutions for a Non-zero Bow Plane 

Angle = - 20 degrees) 
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Figure 31. Degree of Stability for a Non-zero Bow Plane Angle 

= - 20 degrees) 

B. SIMULATIONS USING NUMERICAL INTEGRATION METHODS 

The linearized dynamic response results were verified by simulations using numerical 
integration of the full six degrees of freedom equations of motion for the swimmer delivery 
vehicle (SDV). Figure 32 shows a plot of angle of pitch (0) versus time for the center of 
gravity forward of the center of buoyanc\' case (x^g = + 1) wiih a dive plane angle (6^.) of 

- 15 degrees. The dotted line shows the linearized results from figure 7, the solid line 
shows the numerical integration results. The steady state results of the numerical 
integration method match the linearized dynamic results exactly. 
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Figure 32. Numerical Integration Solution for Angle of Pitch (9) when 
Center of Gravity is Forward = + 1%) 

Figure 33 shows a plot of angle of pitch (0) versus time for the center of gravity aft 
of the center of buoyancy case (x^g = - 1) with a dive plane angle (6^) of - 15 degrees. 



numerical integration results. And once again, the steady state results of the numerical 
integration method match the linearized dynamic results exactly. However, this linearized 
dynamic result was for the vertical plane only, the horizontal plane stability analysis 
indicated that this would be an unstable solution (figure 9). The reason for this 
disagreement in the results is investigated by adding an initial angle of roll to the numerical 
integration analysis. Adding a small angle of initial roll (9 q = 1 degree) caused the vehicle 

to steady out at 137 degrees vice 159 degrees as shown in figure 34. This initial roll angle 
also caused a steady state roll angle of 17 degrees as shown in figure 35. In turn, this 
steady state roll angle caused the steady state yaw velocity (r) shown in figure 36 (i.e. 



Again the dotted line shows the linearized results from figure 7, the solid line shows the 
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molion is no longer resiricied to the vertical plane). Figure 37 shows a plot of z versus 
and indicates that the vehicle is taking a helical ascent as dicussed by Booth [Ref 2: 30“^ 
305]. Therefore, the numerical integration solution resulting in a steady state pitch angle c 
159 degrees (figure 33) is unstable in the horizontal plane as predicted by the linearize 
dynamic response analysis. 




Figure 33. .Numerical Integration Solution for Angle of Pitch (0) when 
Center of Gravity is Aft = - 1%) 
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Figure 34 . Numerical Integration Solution for Angle of Pitch (0) when 
= - 19i: and Initial Roll Angle (<f»y) is 1 degree 
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time 

Figure 35. Numerical Integration Solution for Angle of Roll ( 9 ) when 
^GB - ' and Initial Roll .Angle (<p^) is 1 degree 
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Figure 36. Numerical Integration Solution for Yaw Velocity (r) when 

Initial Roll Angle (<p^) is 1 degree 



6000 




Figure 37. Numerical Integration Solution for Heave Velocity (z) when 
Xqjj = • 1% and Initial Roll Angle ((p^) is 1 degree 
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Some question still remains with regards to the measure of stabilty of the 'inverted 
pendulum' solutions predicted by both the linearized dynamic response analysis and the 
numerical integration analysis. The linearized dynamic response analysis predicts a stable 
solution for the case w'hen center of gravity is placed aft of center of buoyancy 
(Xgb = - 1 %) and the dive plane angle (6^) is - 7 degrees (figure 10). The corresponding 

steady state value for pitch angle was 118 degrees (figure 7). A random persistent roll 
disturbance ((p^) was added to the numerical integration model and the results arc shown in 

figure 38. The solid line indicates the results when a small disturbance is added ((p^ 

centered about 0. 1 degrees), the dashed line indicates the results when a large disturbance 
is added ((p^ centered about 1.0 degrees). As expected, the large disturbance caused the 

vehicle to roll over as shown by the resulting angles of roll (<p) in figure 39. However, the 
vehicle continued to remain stable during small constant random disturbances in the 
inverted position as shown by the resulting angles of roll (<p) in figure 40. This indicates 
that indeed these 'inverted pendulum' solutions have a significant measure of stability. 



69 



I«0 ^ llieln 




38. Numerical Integration Solution for Angle of Pitch (0) when 

Xgb = - and a Persistent Roll Disturbance is Added 
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Figure 39. Numerical Integration Solution for Angle of Roll (<p) when 
Xgb = - and a Large Persistent Roll Disturbance is Added 
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Figure 40 . Numerical Integration Solution for .Angle of Roll (<p) when 
Xgy = - l^c and a Small Persistent Roll Disturbance is .Added 
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VI. CONCLUSIONS AND RECOMMENDATIONS. 



The steady state analysis resulted in four possible solutions provided that vehicle 
motion was restricted to the vertical plane. Analyzing the dynamic stability using the 
steady state results as nominal points generally indicates that only one (if any) of the 
four possible solutions will be stable. There are a few cases where two solutions are 
stable, but these cases (certain non-zero bow plane deflection angles) appear to be the 
exception and not the rule. 

The dynamic stability characteristics of submersiblescan be separated with respeet to 
vertical plane motions (u,w,q,0) and horizontal plane motions (v,p,r,<I)). 

It is possible for submersibles to be dynamically stable with respect to vertical plane 
motion in the inverted (belly up) position during ascents ('Inverted Pendulum' 
stabilization). 

'Inverted Pendulum' stabilization is also possible in the horizontal plane. 

Submersibles are able to maintain this inverted orientation (i.e. ascend belly up 
without rolling over) even under some persistent roll excitation. 

As a recommendation, the dynamic stability analysis should be expanded to include 
the case where angle of roll is 180 degrees, and the ease where angle of roll is non- 
zero (i.e. <|) neither equals zero nor 180 degrees). Analyzing the (p =180 degrees cases 
will only involve changing a few signs with regards to trigonometric functions; 
however, analyzing the non-zero cases will require significant effort. 

Furthermore, identifying and characterizing different stability regions over a range of 
variations of the system parameters should be the matter of future research. 
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APPENDIX A 



SIX DEGREE OF FREEDOM EQUATIONS OF MOTION 

Source: Smith, Crane, and Summey [Reference 1:11-16] 

I. LONGITUDINAL (SURGE) EQUATION OF MOTION 

mu- vr + \vq - xq (q^ -f r2) + xq (pq - r) + zq (pr -f q) = 

^pp P“ ^qq q“ + ^rr ^pr 
+ Xu u + Xvvq vvq + Xyp vp -i- Xyr vr 
+ uq ( Xq5s 65 + Xq 5 b 65) + Xr5r ur6f 

+ Xvv v2 -I- Xy,r\y w^- -I- Xv 6 r uv 6 j- - 1 - uvv ( Xi^-5s 65 + Xw'5b ^b). 
+ u2 ( X5s6s + X5b6b ^b^ + - B) sin 6 

+ U~ Aprop 



2 . 



LATERAL (SWAY) EQUATION OF MOTION 
m V + ur - wp + XQ (pq + r) - vq (p2 +r2) + zq (qr - p) = 
Ypp + Yp r -f Ypq pq-f Yqp qr 
+ Yy V + Ypup + Yp ur + Yyq vq + Yy^p wp + Yyyp wr 
-I- Yy UV -I- Yyy,' vw + Ybp 6p 



'^nose 



^tail 

-I- (W-B) cos 0 sincj) 



Cdv h(x) (v-hxr)2 +Cdz b(x) (w-xq)2 dx 

Ucf(x) 
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3. NORMAL (HE.AVE) EQUATION OF MOTION 

m ’ w - uq + vp + XQ (pr - q) + yc (qr + p) - zq (p- + q^F = 
Zqq + Zpp p- + Zpj pr + Zyy r^j 
+ [Zw w + Zquq + Zyp vp + Zyr vr ] 

+ I Zy/ uw + Zyv (^6s ^s'*'^6b ^b)' 

Coy h(x) (v+xr)2 + 00 ^ b(x) (w-xq)2 j dx 

Ucf(x) 



'^nose 



^tail 



+ (W-B) cos 0 cos (j) 



4. ROLL EQUATION OF MOTION 

Ix P + {h -ly) qr + Ixy (pr - q) - lyz (q^ - r^) 

- Ixz (pq + r) + m jQ (w - uq + vp) - zq (v + ur - wp)_ = 
'Kp p +Kr r + Kpq pq + Kqr qr 

r 

+ Ky V + Kp up + Kr ur + Kyq vq + K^p wp + Ky/j- wr 
+ Ky uv + Ky^y v\v^ + {)’ Q W - yg B) cos 0 cos (p 
- (zqW - zgB) cos 0 sin (j) + u2 Kprop 
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5. PITCH EQUATION OF MOTION 

^ + (^x ■ Jy) pr - Ixy + p) + Jyz (p<^ ‘ 

+ Ixz (p“ - r-) - m XQ (w - uq + vp) - zq (u - vr + wp) = 

:Mq q + Mpp P* + Mpr pr + Mrr r-j 

+ Mw'W + Mq uq + Myp vp + Myr vr j 

+ M\v u\v + Myv V- + u- (M5 s 65 + M55 65) 

■^nose 



^tail 



Coy h(x) (v+xr)2 +Cdz b(x) (w-xq)2 ; x dx 

' Ucftx) 



- (xq W - xg B) cos 6 cos q - [zq W - zg B) sin 0 



6. YAW EQUATION OF MOTION 

^z r + (ly - Ix) Pq ■ ^xy (p" ■ q~) ■ ^fyz (p^ + q) 

+ hz (qr - p) - m XQ (v + ur - wp) - yQ (u - vr + wq) = 

Np p +Nr r + Npq pq + Nqr qr 
+ NyV + Np up + Nj- ur + Nyq vq +N\yp wp + Ny^.j- wr 

+ Ny uv + Nyvv vw + N 5 i- u- dj- 

Cgy h(x) (v+xr)- + Cdz b(x) (w-xq)2 j x dx 

Uct^x) 

+ (xQ W - xg B) cos 0 sin (f) + (yQ W - yg B) sin 0 + u- Nprop 



'^nose 



^tail 
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APPENDIX B 



ROTATION SEQUENCE AND EULER ANGLE RATES 

1 . ROTATION SEQUENCE FOR ([), 0 AND ^ 

Smith, Crane, and Summey [Reference 1:8] descibe the transition from body fixed to 
inertial reference frames as follows: 

Since the equations of motion are referred to an axis system that is fixed for the 
SDV (swimmer delivery vehicle), and thus translates and rotates with it, the 
orientation and position of the moving body axis system relative to a fixed inertial 
reference system must be specified. The orientation of the body axis system with 
respect to the inertial reference system is defined by the standard Euler angles 
(yaw), 6 (pitch), and ([) (roll). The rotation sequence from the inertial reference 

system to the body system is ^ ,0 , and (p as shown in Figure B1 taken from Smith, 
Crane, and Summey [Reference 1:18]. 



2. EULER ANGLE RATES FOR AND ^ 

The Euler angle rates used along with the six equations of motion (Appendix A) in 
order to completely determine the motion of the submersible were specified by Smith, 
Crane, and Summey [Reference 1 :20] to be: 

<p = p + q sin (p tan 0 + r cos (p tan 0 

0 = q cos <p - r sin <p 

sin (p cos (p 

= q ^ ^ 

cos 0 cos 0 
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( 1 ) Veh i c le-Cenrered 
Gravi ty-Di rected System 
pa ra 1 I e 1 to inertial 
axis system. 



(2) Hor i zonta I Flight 
Reference System derived 
from XoYqZo by rotation a- 
bout Z through yaw angle 




>om X“Y“2'‘ by rotation 



about Y“ 
ang I e 9 . 



through p i tch 



erence System derived from 
X‘Y'Z' by rotation about X* 
through rol I angle d- 



Figure Bl. Unit Sphere Development of Euler Angles 
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APPENDIX C 



STEADY STATE COMPUTER PROGRAM 



PROGRAM STEADY 

STEADY STATE SOLUTIONS IN THE VERTICAL PLANE 
DIVE PLANE VARIATION 

REAL L, MASS, IX, IY,IZ,IXZ,IYZ,IXY, LAMBDA 

REAL KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , KVQ , KWP , KWR , KV 

REAL KVW,KPN,KDB 

REAL MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , MW , MW , MDS 
REAL MDB,NDRB 

REAL NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , NV 
REAL NVW,NDRS 

DIMENSION X( 9 ) ,BR( 9 ) ,HH( 9 ) ,VECl(9) 

GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 



C 



PI 

WEIGHT^ 

L 

RHO 

G 

CDO 

MASS 

CDZ 

XWW 

XWDS 

XWDB 

XDSDS 

XDBDB 

CDO 

ZW 

ZDS 

ZDB 

MW 

MDS 

MDB 



= 4 . 0*ATAN( 1 . 0 ) 
= 12000.0 
= 17.425 

= 1.94 

= 32.2 

= 0.0057 

=WEIGHT/G 
=0.5*0.5*RHO 



= 1 
= 4 
= 9 
= -l 
= -8 

=-3 
= -2 
=-2 
= 9 
=-l 
= 1 



710E-01*0 

600E-02*0 

660e-03*0 

160E-02*0 

070E-03*0 

CD0*0 

020E-01*0 

270E-02*0 

270E-02*0 

860E-02*0 

113E-02*0 

113E-02*0 



. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 
. 5*RHO* 



L**2 

L**2 

L**2 

L**2 

L**2 

L**2 

L**2 

L**2 

L**2 

L**3 

L**3 

L**3 



OPEN (21, NAME= ' STl . RES ' , STATUS= ' NEW ' ) 
0PEN(22,NAME='ST2.RES' , STATUS= ' NEW ' ) 
0PEN(23,NAME='ST3.RES' , STATUS='NEW' ) 
OPEN ( 24 ,NAME='ST4 .RES' , STATUS= ' NEW ' ) 
0PEN( 31 ,NAME=' COEF.DAT' , STATUS= ' NEW ' ) 
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DEFINE THE LENGTH X, BREADTH BR, AND HEIGHT HH TERMS 

X(l) = -105.9/12.0 
X(2) = -99.3/12.0 
X(3) = -87.3/12.0 
X(4) = -66.3/12.0 
X(5) = 72.7/12.0 
X(6) = 83.2/12.0 
X(7) = 91.2/12.0 
X(8) = 99.2/12.0 
X(9) = 103.2/12.0 

BR(1) = 0.00/12.0 
BR(2) = 8.24/12.0 
BR(3) = 19.76/12.0 
BR(4) = 29.36/12.0 
BR(5) = 31.85/12.0 
BR(6) = 27.84/12.0 
BR(7) = 21.44/12.0 
BR(8) = 12.00/12.0 
BR(9) = 0.00/12.0 

COMPUTE AREA AND CENTROID 



C 



c 



CALL TRAP( 9 ,BR,X,AREA) 

DO 9 1=1,9 

VECl ( I )=X( I ) *BR( I ) 

9 CONTINUE 

CALL TRAP( 9 ,VEC1 ,X,XAA) 

XA=XAA/AREA 

WRITE (*,1002) 

READ (*,*) DSMIND, DSHAXD, IDS 
DSMIN=DSMIND*PI/180 
DSMAX=DSMAXD*P 1/180 
WRITE (*,1001) 

READ (*,*) RATIO 

WRITE (*,1003) 

READ (*,*) DELB 

DELB=DELB*WEIGHT/100 . 0 
WRITE (*,1004) 

READ (*,*) XGB 

XGB=XGB*L/100 . 0 
WRITE (*,1005) 

READ (*,*) ZGB 

WRITE (*,1006) 

READ ( * , * ) XB 

XB=XB*L/100 . 0 
WRITE (*,1007) 

READ (*,*) ZB 

WRITE (31,*) RATIO, DELB, XGB, ZGB, XB , ZB 
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DO 1 1=1, IDS 

DS=DSMIN+(DSMAX-DSMIN) * (I-1)/(IDS-1) 

IF (DELB.EQ.0.0 ) DELB=0 . 0 0 0 0 0 1 
IF (ZGB.EQ.O.O) ZGB =0.000001 
DB=RATIO*DS 

PX =XGB*WEIGHT-XB*DELB 
PZ =ZGB*WEIGHT-ZB*DELB 
DEN =CDZ*AREA* ( PX+XA*DELB ) 

LAMBDA=MW*DELB-PX*ZW+PZ* (XWDS+RATIO*XWDB) *DS 
ALPHA =-PX*( ZDS+RATIO*ZDB) *DS-PZ*CDO+PZ*(XDSDS+ 
RATI O* RATI 0*XDBDB ) *DS**2+DELB* ( MDS+RAT 10 
*MDB) *DS 
BETA =PZ*XWW 
LAMBDA=LAMBDA/DEN 
ALPHA =ALPHA /DEN 
BETA =BETA /DEN 
C 
C 

A = 1.0+BETA 

B = LAMBDA 

C = ALPHA 

DET= B**2-4.0*A*C 

IF (DET.LT.0.0) GO TO 2 

WP=( -B+SQRT( DET) )/(2.0*A) 

YY=-XWW*WP * * 2- ( XWDS*DS+XWDB*DB ) *WP 
& -(XDSDS*DS**2+XDBDB*DB**2)+CD0 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
& *WP*ABS(WP) 

IF (WP.LT.0,0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
& *WP*ABS(WP) 

THETA=ATAN2 ( YY , XX ) 

USQ=DELB * SIN ( THETA )/YY 
THETA=THETA*180/PI 
DSD=DS*180/PI 
IF (USQ.LT.0.0) GO TO 3 
IF (WP.GE.0,0) U= SQRT(USQ) 

IF (WP.LT.O.O) U=-SQRT(USQ) 

W=WP*U 

WRITE (21,*) DSD, THETA, U,W,WP 
C 

3 WP=(-B-SQRT(DET) )/(2.0*A) 

YY=-XWW*WP**2-(XWDS*DS+XWDB*DB ) *WP 
& -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF (WP.GE.0,0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
& *WP*ABS(WP) 

IF (WP.LT.O.O) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
& *WP*ABS(WP) 

THETA=ATAN2 ( YY , XX ) 

USQ=DELB * SIN ( THETA )/YY 
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DSD=DS*180/PI 
THETA=THETA*180/PI 
IF (USQ.LT. 0.0) GO TO 2 
IF (WP.LT.0.0) U=-SQRT(USQ) 

IF (WP.GE.0.0) U= SQRT(USQ) 
W=WP*U 

WRITE (22,*) DSD, THETA, U,W,WP 



2 A = -1.0+BETA 

DET= B**2-4.0*A*C 
IF (DET.LT.0.0) GO TO 1 
WP=(-B+SQRT(DET) )/(2. 0*A) 
YY=-XWW*WP**2-(XWDS*DS+XWDB*DB ) *WP 
& -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS(WP) 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS (WP) 

THETA=ATAN2 ( YY , XX ) 

USQ=DELB*S IN ( THETA )/YY 
DSD=DS*180/PI 
THETA=THETA*180/PI 
IF (USQ.LT. 0.0) GO TO A 
IF (WP.GE.0.0) U— SQRT(USQ) 

IF (WP.LT.0.0) U= SQRT(USQ) 

W=WP*U 

WRITE (23,*) DSD, THETA, U,W,WP 
4 WP=(-B-SQRT(DET) )/(2 . 0*A) 

YY=-XWW*WP**2-( XWDS*DS+XWDB*DB ) *WP 
& - ( XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS ( WP ) 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS (WP) 

THETA=ATAN2 ( YY , XX ) 

USQ=DELB *S IN ( THETA )/YY 
DSD=DS*180/PI 
THETA=THETA*180/PI 
IF (USQ.LT. O.O) GO TO 1 
IF (WP.GE.0.0) U— SQRT(USQ) 

IF (WP.LT.0.0) U* SQRT(USQ) 

W=WP*U 

WRITE (24,*) DSD, THETA, U,W,WP 



1 


CONTINUE 








STOP 








1001 


FORMAT 


( ' 


ENTER 


BOW PLANE TO DIVE PLANE RATIO' 


1002 


FORMAT 


( ' 


ENTER 


MIN, MAX, AND INCREMENTS IN 


& 






DS ( degrees ) ' ) 


1003 


FORMAT 


( ' 


ENTER 


DELB ( %W) ' ) 


1004 


FORMAT 


( ' 


ENTER 


XGB ( %L ) ' ) 



SI 



noon 



1005 FORMAT (' ENTER ZGB (feet)') 

1006 FORMAT (' ENTER XB (%L)') 

1007 FORMAT (' ENTER ZB (feet)') 

END 

C 

SUBROUTINE TRAP ( N , A , B , OUT ) 

NUMERICAL INTEGRATION ROUTINE USING 
THE TRAPEZOIDAL RULE 

DIMENSION A( 1 ) ,B( 1 ) 

N1=N-1 
OUT=0 . 0 
DO 1 1=1, Nl 

OUT1=0 . 5* ( A( I )+A( I+l ) )*(B(I+1)-B(I) ) 
OUT =OUT+OUTl 
1 CONTINUE 
RETURN 
END 
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APPENDIX D 



LINEARIZED DYNAMIC STABILITY COMPUTER PROGRAM 



PROGRAM LINEARIZED DYNAMIC STABILITY 

10 20 30 40 50 

345678901234567890123456789012345678901234567890123456 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L , MASS , lA, IX, I Y, I Z 

DOUBLE PRECISION KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , 

& KVQ,KWP,KWR,KV, 

& KVW , KPN , KDB , MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , 

& MW,MW,MDS,MDB, 

& NDRB , NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , 

& NV,NVW,NDRS 

DIMENSION A1(4,4),B1(4,4) ,BETAl ( 4 ) ,ALFR1 ( 4 ) ,ALFI1 ( 4 ) 
DIMENSION BB1(4,4),BB2(4,4),ZZZ1(4,4),ZZZ2(4,4) 
DIMENSION A2(4,4),B2(4,4), BETA2 ( 4 ) ,ALFR2 ( 4 ) , ALFI2 ( 4 ) 
DIMENSION WRl ( 4 ) ,WR2 ( 4 ) , KI 1 ( 4 ) , WI 2 ( 4 ) 

DIMENSION X(9) ,BR(9) ,VEC1(9) ,VEC2(9) 

GEOMETRIC PROPERTIES 



PI = 4 .D0*DATAN(1 .DO) 

WEIGHT= 12000.0 

IX = 1760.0 

lY =9450.0 

IZ =10700.0 

L = 17.425 

RHO = 1.94 

G = 32.2 

CDO = 0.0057 

MASS = KEIGHT/G 

CDZ = 0.5*0.5*RHO 

CDO = CD0*0. 5*RHO*L**2 

SURGE HYDRODYNAMIC COEFFICIENTS 

XPP = 7 . 030E-03*0 . 5*RHO’^L**4 

XQO =-l . 470E-02*0 . 5’^RHO*L** 4 

XRR = 4 . 010E-03*0 . 5*RHO*L**4 

XPR = 7 .640E-04*0. 5*RHO*L**4 

XUDOT =-7 . 580E-03*0 . 5*RHO*L* *3 

XWQ =-l . 920E-01*0 . 5*RH0*L**3 
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c 

c 

c 



c 

c 

c 



c 

c 

c 



XVP =-3 
XVR = 1 
XQDS = 2 
XQDB =-2 
XRDR =-8 
XW = 5 
XWW = 1 
XVDR = 1 
XWDS = 4 
XWDB = 9 
XDSDS =-l 
XDBDB =-8 
XDRDR =-l 
XRES = 



, 240E-03*0 
.890E-02*0 
. 610E-02*0 
.600E-03*0 
. 180E-04*0 
.290E-02*0 
. 710E-01*0 
.730E-03*0 
.600E-02*0 
. 660E-03*0 
.160E-02*0 
. 070E-03*0 
. 010E-02*0 
CD0*0 



5*RHO*L*’"3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 



5*RHO*L**2 



SWAY HYDRODYNAMIC COEFFICIENTS 



YPDOT = 1 
YRDOT = 1 
YPQ = 4 
YQR =-6 
YVDOT =-5 
YP =3 
YR =2 
YVQ = 2 
YWP = 2 
YWR =-l 
YV =-9 
YW = 6 
YDRS =+2 
YDRB =+2 



270E-04*0 . 
240E-03*0 . 
125E-03*0 . 
510E-03*0. 
550E-02*0 . 
055E-03*0 . 
970E-02*0 . 
360E-02*0 . 
350E-01*0 . 
880E-02*0 . 
310E-02*0 . 
840E-02*0 . 
270E-02*0 . 
270E-02*0 . 



5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 



HEAVE HYDRODYNAMIC COEFFICIENTS 



ZQDOT 


=-6.810E- 


ZPP 


= 1.270E- 


ZPR 


= 6.670E- 


ZRR 


=-7 . 350E- 


ZWDOT 


=-2.430E- 


ZQ 


=-l . 350E- 


ZVP 


=-4.810E- 


ZVR 


= 4.550E- 


zw 


=-3 . 020E- 


zvv 


=-6.840E- 


ZDS 


=-2 . 270E- 


ZDB 


=-2.270E- 



03*0 . 5*RHO*L**4 
04*0 . 5*RHO*L**4 
03*0 . 5*RHO*L**4 
03*0 . 5*RHO*L**4 
01*0 . 5*RHO*L**3 
01*0 . 5*RHO*L**3 
02*0 . 5*RHO*L**3 
02*0 . 5*RHO*L**3 
01*0 . 5*RHO*L**2 
02*0 . 5*RHO*L**2 
02*0 . 5*RHO*L**2 
02*0 . 5*RHO*L**2 



ROLL HYDRODYNAMIC COEFFICIENTS 



KPDOT =-l . 010E-03*0 . 5*RH0*L**5 
KRDOT =-3 . 370E-05*0 . 5*RHO*L**5 
KPQ =-6 . 930E-05*0 . 5*RHO*L**5 
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KQR 


= 


1 . 680E-02*0 . 5*RHO*L**5 


KVDOT 


= 


1.270E-04*0. 5*RHO*L**4 


KP 


= 


-1. 100 E- 02*0. 5*RHO*L**4 


KR 


= 


-8.410E-04*0. 5*RHO*L**4 


KVO 




-5 . 115E-0 3*0 . 5*RHO*L**4 


KWP 


= 


-1 . 270E-04*0 . 5*RHO*L**4 


KWR 


= 


1 . 3 90E-02*0 . 5*RHO*L**4 


KV 


= 


3 . 0 55E-03*0 . 5*RHO*L**3 


KVW 


= 


-1. 870 E- 01*0. 5*RH0*L**3 


PITCH 


HYDRODYNAMIC COEFFICIENT 


MQDOT 




-1.680E-02*0. 5*RHO*L**5 


KPP 




5 . 260E-05*0 . 5*RH0*L**5 


NPR 


= 


5.040E-03*0. 5*RHO*L**5 


MRR 




-2.860E-03*0. 5*RHO*L**5 


MWDOT 


= 


-6 . 810E-02*0 . 5*RHO*L**4 


MQ 


= 


-6.860E-02*0. 5*RHO*L**4 


MVP 


= 


1.180E-03*0. 5*RH0*L**4 


MVR 


= 


1. 730 E- 02*0. 5*RHO*L**4 


MW 


= 


9.860E-02*0. 5*RH0*L**3 


MW 


= 


-2, 510 E- 02*0. 5*RHO*L**3 


MDS 


= 


-1. 113 E- 02*0. 5*RHO*L**3 


MDB 




1. 113 E- 02*0. 5*RHO*L**3 


YAW HYDRODYNAMIC COEFFICIENTS 


NPDOT 




-3 . 37 0E-05*0 . 5*RHO*L**5 


NRDOT 


= 


-3. 400 E- 03*0. 5*RHO*L**5 


NPQ 


= 


-2.110E-02*0 . 5*RHO*L**5 


NQR 


= 


2.750E-03*0. 5*RHO*L**5 


NVDOT 


= 


1 . 24 0E-03*0 . 5*RHO*L**4 


NP 


- 


-8. 405 E- 04*0. 5*RHO*L**4 


NR 


= 


-1. 640 E- 02*0. 5*RHO*L**4 


NVQ 


= 


-9 . 990E-03*0 . 5*RHO*L**4 


NWP 


= 


-1 . 750E-02*0 . 5*RHO*L**4 


NWR 


= 


7. 350 E- 03*0. 5*RHO*L**4 


NV 


= 


-7. 420 E- 03*0. 5*RHO*L**3 


NVW 


= 


-2. 670 E- 02*0. 5*RHO*L**3 


NDRS 


= 


-1. 113 E- 02*0. 5*RHO*L**3 


NDRB 


= 


+1. 113 E- 02*0. 5*RHO*L**3 


DEFINE 


THE LENGTH X AND BREADTH BR TERMS 


X(1 ) 




-105.9/12.0 


X(2) 


= 


-99.3/12.0 


X( 3 ) 


= 


-87 . 3/12.0 


X( 4 ) 


= 


-66.3/12.0 


X(5) 


= 


72.7/12.0 


X(6) 


= 


83.2/12.0 


X(7) 


= 


91.2/12.0 


X(8) 


= 


99.2/12.0 


X(9) 


= 


103 . 2/12 . 0 
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c 



BR{1) 
BR( 2 ) 
BR(3) 
BR( 4 ) 
BR( 5) 
BR( 6 ) 
BR(7) 
BR(8) 
BR( 9 ) 



0 . 00 / 12.0 

8.24/12.0 

19.76/12.0 

29.36/12.0 

31.85/12.0 

27.84/12.0 

21.44/12.0 

12 . 00 / 12.0 

0 . 00 / 12.0 



COMPUTE AREA, CENTROID, AND MOMENT OF INERTIA 

CALL TRAP( 9 ,BR,X,AREA) 

DO 9 1=1,9 

VECl ( I )=X( I )*BR( I ) 

VEC2(I)=X(I)*VEC1(I) 

9 CONTINUE 

CALL TRAP( 9 ,VEC1 ,X,XAA) 

XA=XAA/AREA 

CALL TRAP( 9 , VEC2 ,X, lA) 



WRITE (*,1001) 

READ (*,*) IRES 

OPEN( 31 ,NAME= ' COEF.DAT' , STATUS= 'OLD' ) 

READ(31,*) RATIO, DELB, XGB, ZGB, XB , ZB 

BUOY= WEIGHT + DELB 

XG=XB+XGB 

ZG=ZB+ZGB 



MASS MATRIX COEFFICIENTS 



C 



C 



C 



C 



Bl(l,l)= MASS - XUDOT 
Bl(l,2)= 0.0 
Bl ( 1 , 3 )= MASS*ZG 
Bl(l,4 )= 0.0 

Bl(2,l)= 0.0 
Bl(2,2)= MASS - ZWDOT 
Bl ( 2 , 3 ) =- ( ZQDOT+MASS*XG ) 
Bl ( 2 , 4 )= 0.0 

Bl ( 3 , 1 )= MASS*ZG 
Bl( 3,2)=-(MWDOT+MASS*XG) 
Bl(3,3)= lY-MQDOT 
Bl(3,4)= 0.0 



Bl( 4 ,1)= 0.0 
Bl(4,2)= 0.0 
Bl(4,3)= 0.0 
Bl ( 4 , 4 )= 1.0 
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B2(l,l)= IX-KPDOT 
B2(l,2)= 0.0 

B2 ( 1 , 3 )=- ( KVDOT+MASS*ZG) 

B2 ( 1 , <3 )=-KRDOT 
C 

B2(2,1 )= 0.0 
B2(2,2)= 1.0 
B2(2,3)= 0.0 
B2(2,4 )= 0.0 
C 

B2(3,1 )=-(YPD0T+KASS*ZG) 

B2(3,2)= 0.0 
B2(3,3)= MASS-YVDOT 
B2(3,4)= MASS*XG-YRDOT 
C 

B2 ( 4 , 1 )=-NPDOT 
B2(4,2)= 0.0 
B2(4,3)= MASS*XG-NVD0T 
B2(4,4)= IZ-NRDOT 
C 

OPEN( 41 ,NAME= ' DECS . RES ' , £TATUS= 'NEW' ) 

OPEN( 42,NAME='DE0S1 .RES' , STATUS='NEK' ) 
OPEN(43,NAKE='DEOS2.RES' , STATUS= ' NEW ' ) 

C 

IF ( IRES . EQ. 1 ) GO TO 1 
IF ( IRES.EQ.2 ) GO TO 2 
IF (IRES.EQ.3) GO TO 3 
IF ( IRES.EQ.4 ) GO TO 4 

1 OPEN ( 21 ,NAME= ' STl . RES' , STATUS= 'OLD' ) 

11 READ ( 21 , * , END=100 ) DSD , THETO , UO , WO , WP 
GO TO 5 

2 OPEN ( 22,NANE='ST2 .RES' ,STATUS='OLD' ) 

12 READ ( 22 , * , END=100 ) DSD , TKETO , UO ,W0 , WP 
GO TO 5 

3 OPEN (23,NAME='ST3.RES' ,STATUS='OLD' ) 

13 READ ( 23 , * , END=100 ) DSD , THETO , UO , WO , WP 
GO TO 5 

4 OPEN ( 24,NAKE='ST4 .RES' ,STATUS='OLD' ) 

14 READ ( 24 , * , END=100 ) DSD , THETO , UO , WO , WP 
GO TO 5 

C 

5 THETA0=THET0*PI/180 . 0 
DS = DSD*PI/180.0 

DB = DS * RATIO 
C 

C DAMPING MATRIX COEFFICIENTS 

C 

A1 ( 1 , 1 )=-2 . 0*U0*CD0+W0^ ( XWDS-DS-^XWDB''-DB ) 

& +2.0*U0* (XDSDS^DS’"*2-XD5DB*DE**2 ) 

A1 ( 1 , 2 )= 2 . 0*XWTn’*W 0 + U0* ( XWDS-DS + XWDB’^DB ) 
A1 { 1 , 3 )= (XWQ-MiASS )*W0^ (XODS-DS-^XQDB’^DB) ’^UO 
Al ( 1 , 4 )=- (WEIGHT-BUOY ) -DCOS ( THETAO ) 
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c 



c 



c 



c 



c 



c 



c 



A1 ( 2 , 1 )= ZW*W0+2 . 0*U0* ( ZDS*DS+ZDB*DB ) 

A1 ( 2 , 2 )- ZW*U0-2 . 0*CDZ*AREA*DABS(W0 ) 

A1 ( 2 , 3 )= ( ZQ+MASS ) *U0+2 . 0*CDZ*AREA*XA*DABS ( WO ) 
Al ( 2 , 4 )=- (WEIGHT-BUOY) *DS IN ( THETAO ) 



& 

Sc 



Al(3,l)= 

Al(3,2)= 

Al(3,3)= 

Al(3,4)= 



MW*W0+2 . 0*U0* (MDS*DS+MDB*DB ) 
MW*U0+2 . 0*CDZ*AREA*XA*DABS(W0 ) 
(MQ-MASS*XG) *U0-MASS*ZG*W0 
-2. 0*CDZ*IA*DABS(W0 ) 

( XG*WEIGHT-XB*BUOY ) *DS IN ( THETAO ) - 
( ZG*WEIGHT-ZB*BUOY ) *DCOS ( THETAO ) 



Al(4,l)= 0.0 
Al(4,2)= 0.0 
Al( 4 ,3)= 1.0 
Al(4,4)= 0.0 



A2(l,l)= KP*UO+(KWP-MASS*ZG)*WO 

A2 ( 1 , 2 )=- ( ZG*WE I GHT-ZB *BUOY ) *DCOS( THETAO ) 

A2(l,3)= KV*U0+KVW*W0 

A2(l,4)= (KR+MASS*ZG)*UO+KWR*WO 



A2(2,l)= 
A2 ( 2 , 2 ) = 
A2(2,3)= 
A2 ( 2 , 4 ) = 



1 . 0 
0 . 0 
0 . 0 

DTAN( THETAO ) 



A2(3,l)= 
A2(3,2)= 
A2( 3 , 3 ) = 
A2(3,4)= 



YP*U0+( YWP+MASS ) *W0 
( WEIGHT-BUOY ) *DCOS ( THETAO ) 
YV*UO+YVW*WO-CDZ*AREA*DABS(WO ) 
YWR*W0+ ( YR-KASS ) *UO-CDZ*AREA*XA 
*DABS(W0 ) 



A2(4,l)= 
A2 ( 4 , 2 ) = 
A2(4,3)= 
A2(4,4)= 



MASS*XG*W0+NP*U0+NWP*W0 
( XG*WEIGHT-XB*BUOY ) *DCOS ( THETAO ) 
NV*UO+NVW*WO-CDZ*AREA*XA*DABS (WO ) 
(NR-MASS*XG) *UO+NWR*WO-CDZ*IA*DABS (WO ) 



RESTORE B-MATRIX 



DO 71 1=1,4 
DO 72 J=l,4 

BB1(I,J)=B1(I,J) 
72 CONTINUE 

71 CONTINUE 

DO 81 1=1,4 
DO 82 J=l,4 

BB2(I,J)=B2(I,J) 
82 CONTINUE 

81 CONTINUE 
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c 

CALL RGG( 4 , 4 , Al , BBl , ALFRl ,ALFIl , BETAl , 0,ZZZl , lER) 
CALL DEGSTB(DE0S1 , ALFRl ,ALFIl , BETAl , FREQl ,WRl ,V7Il ) 
C 

CALL RGG( 4 , 4 , A2 , BB2 , ALFR2 ,ALFI2 , BETA2 , 0 , ZZZ2 , lER ) 
CALL DEGSTB( DE0S2 ,ALFR2 , ALFI2 , BETA2 , FREQ2 ,WR2 ,WI2 ) 
C 

IF (DEOSl .GE.DE0S2 ) DEOS=DEOSl 
IF ( DEOSl . LT.DEOS2 ) DEOS-DEOS2 
C 

WRITE (41,2001) DSD,THETO,UO,WO,WP,DEOS, 

& DEOSl, DEOS2 

IF ( DEOS . LT. 0 . DO ) 

& WRITE (42,2001) DSD,THETO,UO,WO,WP,DEOS, 

& DEOSl, DEOS2 

IF (DEOSl .LT.O. DO) 

& WRITE (43,2001) DSD , THETO , UO , WO , WP , DEOS , 

& DEOSl, DEOS2 

C 

IF ( IRES.EQ. 1 ) GO TO 11 
IF ( IRES.EQ, 2 ) GO TO 12 
IF ( IRES . EQ, 3 ) GO TO 13 
IF ( IRES.EQ. 4 ) GO TO 14 
C 

100 STOP 
1001 FORMAT 
& 

2001 FORMAT 

2002 FORMAT 
END 

C 

SUBROUTINE DEGSTB ( DEOS , ALFR , ALF I , BETA , OMEGA , WR , WI ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION ALFR( 4 ) , ALFI ( 4 ) , BETA( 4 ) , WR ( 4 ) , WI ( 4 ) 

DO 1 1=1,4 

WR( I )=ALFR( I )/BETA( I ) 

WI ( I )=ALFI ( I )/BETA( I ) 

1 CONTINUE 
DEOS=-l . OE+10 
DO 2 1=1,4 

IF ( WR( I ) , LT.DEOS ) GO TO 2 
DEOS=WR( I ) 

U = I 

2 CONTINUE 
OMEGA=KI ( I J ) 

OMEGA=DABS ( OMEGA ) 

RETURN 
END 

C 

SUBROUTINE TRAP ( N , A , B , OUT ) 



( ' ENTER THE RESPONSE DATA FILE DESIRED 
(1,2,3, OR 4) ' ) 

(8E15.5) 

(F10.3) 
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NUMERICAL INTEGRATION ROUTINE USING 
THE TRAPEZOIDAL RULE 

IMPLICIT DOUBLE PRECISION (A-H,0-2) 
DIMENSION A( 1 ) ,B( 1 ) 

Nl=N-l 
OUT=0 . 0 
DO 1 1=1, N1 

OUT1 = 0 . 5*(A(I)+A(I + 1))*(B(I + 1)-B(D) 
OUT =OUT+OUTl 
1 CONTINUE 
RETURN 
END 



90 



LIST OF REFERENCES 



1. Smith, N.S.. Crane, J.W., and Summey,D.C.,"SDV Simulator Hydrodynamic 
Coefficients." NCSC Technical Memorandum 231-78, June 1978. 

2. Booth. T.B.. "Stability of Buoyant Underwater Vehicles, Predominantly Forward 
Motion." International Shipbuilding Progress, Volume 24, November 1977, 
Nr. 279, pages 297-305. 

3. Booth, T.B.. "Stability of Buoyant Underwater Vehieles, Near Vertical Ascent," 
International Shipbuilding Progress, Volume 24, December 1977, Nr. 280, 
pages 346-352. 



BIBLIOGRAPHY 



Boncal, R.J., "A Study of Model Based Manuevering Controls for Autonomous 
Underwater Vehicles," Master's Thesis, Naval Postgraduate School, Monterey, 
California, December 1987 

Gueler. G.F., "Modelling. Design, and Analysis of an Autopilot for Submersible 
Vehicles." International Shipbuilding Progress, Volume 36, February 1989, 
Nr. 405, pages 51-85. 

Richards, R.J. and Stolen, D.P., "Depth Control of a Submersible Vehicle," 
International Shipbuilding Progress, Volume 28, February 1981, Nr. 318, 
pages 30-39. 



91 



INITIAL DISTRIBUTION LIST 



No. Copies 

1 . Defense Technical Information Center 2 

Cameron Station 

Alexandria, Virginia 22304-6145 

2. Library, Code 52 2 

Naval Postgraduate School 

Monterey, California 93943-5002 

3. Department Chairman, Code ME 1 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 93943-5000 

4. Professor Fotis A. Papoulias 2 

Code ME/Pa 

Department of Mechanical Engineering 
Naval Postgraduate School 
Monterey, California 93943-5000 

5. LT Brian D. McKinley 2 

Navy Experimental Diving Unit 

Panama City, Florida 32407 

6. Naval Engineering Curricular Officer, Code 34 1 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 93943-5004 



92 



DUDLEY KNO- -JSRARY 
NAVAL pr^S : GRADUATE SCHOOL 
MONTERt;>‘ CA 93943-5101 



GAYLORD S 



